{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a13b877f-af1d-472e-8e93-b917ccc0a018",
   "metadata": {},
   "outputs": [],
   "source": [
    "import calfem.geometry as cfg\n",
    "import calfem.mesh as cfm\n",
    "import calfem.vis as cfv\n",
    "import numpy as np\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8613b582-bae4-42dc-8b70-9b723c93f84c",
   "metadata": {},
   "source": [
    "# Define the geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a5e2d391-1bed-468f-bd59-a4192b8ce0c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "g = cfg.Geometry()\n",
    "\n",
    "g.point([0.0, 0.0]) # point 0\n",
    "g.point([1.0, 0.0]) # point 1\n",
    "g.point([1.0, 1.0]) # point 2\n",
    "g.point([0.0, 1.0]) # point 3\n",
    "\n",
    "# g.spline([0, 1]) # line 0\n",
    "# g.spline([1, 2]) # line 1\n",
    "# g.spline([2, 3]) # line 2\n",
    "# g.spline([3, 0]) # line 3\n",
    "# g.surface([0, 1, 2, 3])\n",
    "\n",
    "factor = 50\n",
    "dir_id_1 = 100\n",
    "dir_id_2 = 200\n",
    "dir_id_3 = 300\n",
    "dir_id_4 = 400\n",
    "\n",
    "g.spline([0, 1], el_on_curve=factor, marker=dir_id_1) # line 0\n",
    "g.spline([1, 2], el_on_curve=factor, marker=dir_id_2) # line 1\n",
    "g.spline([2, 3], el_on_curve=factor, marker=dir_id_3) # line 2\n",
    "g.spline([3, 0], el_on_curve=factor, marker=dir_id_4) # line 3\n",
    "g.structuredSurface([0, 1, 2, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86e70042-4dce-41ca-99a0-7820bd27d437",
   "metadata": {},
   "source": [
    "# Meshing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4ac3b026-7104-44c2-801a-170e3e71025d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Info    : GMSH -> Python-module\n"
     ]
    }
   ],
   "source": [
    "cfv.drawGeometry(g)\n",
    "cfv.showAndWait()\n",
    "\n",
    "mesh = cfm.GmshMesh(g,return_boundary_elements=True)\n",
    "\n",
    "mesh.elType = 2          # Degrees of freedom per node.\n",
    "mesh.dofsPerNode = 1     # Factor that changes element sizes.\n",
    "# mesh.elSizeFactor = 0.5 # Element size Factor\n",
    "\n",
    "coords, edof, dofs, bdofs, elementmarkers, boundaryElements  = mesh.create()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "19896c03-b142-4397-8910-0ec1d4bc5f64",
   "metadata": {},
   "outputs": [],
   "source": [
    "cfv.figure()\n",
    "\n",
    "# Draw the mesh.\n",
    "\n",
    "cfv.drawMesh(\n",
    "    coords=coords,\n",
    "    edof=edof,\n",
    "    dofs_per_node=mesh.dofsPerNode,\n",
    "    el_type=mesh.elType,\n",
    "    filled=True,\n",
    "    title=\"Example 01\"\n",
    "        )\n",
    "cfv.showAndWait()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26faf75c-b847-42d3-aa53-7faf28661c27",
   "metadata": {},
   "source": [
    "# Generate element, node and edge, normal vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d28594f7-1a70-4771-8239-5a1947b8cc7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "element = edof - 1\n",
    "node = coords\n",
    "\n",
    "E = len(element)\n",
    "N = len(node)\n",
    "\n",
    "dir_edge_1 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_1]])-1\n",
    "dir_edge_2 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_2]])-1\n",
    "dir_edge_3 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_3]])-1\n",
    "dir_edge_4 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_4]])-1\n",
    "\n",
    "# boundary edge with node in counterclock-wise order\n",
    "boundary = [edge['node-number-list'] for marker in boundaryElements for edge in boundaryElements[marker]]\n",
    "\n",
    "edge = set([tuple(x) for x in boundary])\n",
    "\n",
    "# triangle - E x 3\n",
    "eta = np.ones((E,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3df49b0b-9c42-4c48-8939-8932b2a1ad89",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, e in enumerate(element):\n",
    "    if (e[0],e[1]) not in edge:\n",
    "        if (e[1],e[0]) not in edge:\n",
    "            edge.add((e[0],e[1]))\n",
    "        else:\n",
    "            eta[i,0] = -1\n",
    "    if (e[1],e[2]) not in edge:\n",
    "        if (e[2],e[1]) not in edge:\n",
    "            edge.add((e[1],e[2]))\n",
    "        else:\n",
    "            eta[i,1] = -1            \n",
    "    if (e[2],e[0]) not in edge:\n",
    "        if (e[0],e[2]) not in edge:\n",
    "            edge.add((e[2],e[0]))\n",
    "        else:\n",
    "            eta[i,2] = -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "66619742-18e1-4ca3-a14d-63b47c196552",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.  0.  0. ...  0.  0.  0.]\n",
      " [ 0.  1.  0. ...  0.  0.  0.]\n",
      " [ 0.  0.  1. ...  0.  0.  0.]\n",
      " ...\n",
      " [ 0.  0.  0. ...  4. -1.  0.]\n",
      " [ 0.  0.  0. ... -1.  4. -1.]\n",
      " [ 0.  0.  0. ...  0. -1.  4.]]\n",
      "2600\n",
      "(2601, 2601)\n"
     ]
    }
   ],
   "source": [
    "N = len(node)\n",
    "F = np.zeros((N,N))\n",
    "\n",
    "for i, e in enumerate(element):\n",
    "    x1,y1 = node[e][0]\n",
    "    x2,y2 = node[e][1]\n",
    "    x3,y3 = node[e][2]\n",
    "    \n",
    "    A_e = 1/2*np.linalg.det(np.c_[node[e], np.ones((3,1))])\n",
    "    \n",
    "    L1 = np.sqrt((x2-x1)**2+(y2-y1)**2)\n",
    "    L2 = np.sqrt((x3-x2)**2+(y3-y2)**2)\n",
    "    L3 = np.sqrt((x1-x3)**2+(y1-y3)**2)\n",
    "    \n",
    "    left = np.linalg.inv([[(y2-y1)/L1,(x1-x2)/L1],\n",
    "                          [(y3-y2)/L2,(x2-x3)/L2]])\n",
    "    \n",
    "    Phi_L = np.c_[left, np.zeros((2,1))]\n",
    "    Phi_e = Phi_L  @ np.array([[-1/L1,1/L1,0],\n",
    "                               [0,-1/L2,1/L2],\n",
    "                               [1/L3,0,-1/L3]])\n",
    "    \n",
    "    Phi_e = (1/(2*A_e))*np.array([[x2-x3,x3-x1,x1-x2],\n",
    "                                  [y2-y3,y3-y1,y1-y2]])\n",
    "    \n",
    "    # print(np.linalg.norm(Phi_e-_Phi_e))\n",
    "    \n",
    "    F_e = A_e * Phi_e.T @ Phi_e\n",
    "    \n",
    "    F[np.ix_(e,e)] = F[np.ix_(e,e)] + F_e\n",
    "    # with np.printoptions(precision=2, suppress=True):\n",
    "    #     print(e)\n",
    "    #     print(F_e)\n",
    "    #     print(F)\n",
    "    \n",
    "print(F)\n",
    "print(np.linalg.matrix_rank(F))\n",
    "print(F.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "6b52107d-967b-4e26-b7e0-87b8bccf392d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x205431e9f70>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQwAAAD8CAYAAACCaZo+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAATs0lEQVR4nO3dX6wcZ33G8e9zjp0gTFCchlqubZWA3Atz0WAsJyoIUSEcxzcONyi5IFYayVw4FUj0woGLIFAkWhWQItFIRlg4FSWKBCgWchuMRRVxkeADCk6cNPEhfxRbji3qKKSJFP85v17se+w9692zs7uzs+/MPh9ptbPvzO6+Z7zzeN7fznmPIgIzsyJmJt0BM6sPB4aZFebAMLPCHBhmVpgDw8wKc2CYWWGVB4ak7ZJelDQvaW/V729mw1OV12FImgVeAj4HnASOAndFxPOVdcLMhlb1GcZWYD4iXo6I88CjwM6K+2BmQ1pR8futA15ve3wSuKV9A0m7gd0AM7PXfOL97//QlXVA+/nQqI+7qeI9zHLyf2+f+lNEfKj/ltUHRl8RsQ/YB3DdB9fHJ/7uH9HF5hyCiiCkno+HeU7Zj226/PeR+18rum3VQ5JTwIa2x+tTW0+6GMSK5nyYOw/MIgdqv+eU/bgM6qiNlf3YJqPqwDgKbJR0k6RrgDuBg/2e1LTQmAYOtXL6kJtKhyQRcVHSfcATwCywPyKO99q+/SOyGBpNGp7YZOUQaqO+xjDDTwJiZrgArryGERGHgEOFtu147NAwW2qYUJtZWICF4UKjdld6enhiNpqF2RkUgRYG/4+3doEBDg2zUQ0bGlkHxnKR4NAwG60wuxgag8g6MPpeZOXQsBLl8K3JoEYtzC7MDhYBWQdGEQ6NPI374BvHwZjDtya5q31gQL1Co4wDwQdfMw/GOsg6MAb5SBQJDS3EkoOjioOzUxkHgg8+m5SsA2PQ//eWC43FanD7wVHFwWnWJFkHxjC6hcblsBjy6jYza2lcYMDS0HBYmJUn68AY5RDXxSBmWi+yGBaKqOSrLrOmyjowRjm0tRDMXAouvm92SfJ0Fj7NrLisA2NY7QXO2fMLXLpmpnWmIREzcmiYDalxgXFVzSJwaJiVJOvAGLSG0RkWl2sWDg2zUmQdGAMfyrr625DLoeDQMBtZ1oExqG4XUS0JBYeG2UgaFRjdODTMytOowOh1wDs0zMrRqMBY7vc4HBpmo2tUYPTj0DAbzVQFBjg0zEYxdYEBDg2zYWUdGOP8/VKHhtngsg6McR+mDg2zwWQdGFVwaJgVN/WBAQ4Ns6KyDoyBf/lshIPZoWHWX9aBMehhOeoEvA4Ns+VlHRiT4NAw682B0YVDw6y7rAOjyhpGp56hca1Dw6ZX1oFRdQ2j2+tdFRrv+UzDptdIgSHpVUnPSnpG0lxqu0HSYUkn0v3q1C5JD0mal3RM0uYyfoBx8/DE7IoyzjD+PiJujogt6fFe4EhEbASOpMcAtwMb02038HAJ710Jh4ZZyziGJDuBA2n5AHBHW/sj0fIUcL2ktcu90CRrGJ1c0zAbPTAC+KWk30nandrWRMTptPwGsCYtrwNeb3vuydS2hKTdkuYkzZ2/8M5gnRnzH0J2TcOm3YoRn/+piDgl6S+Bw5L+p31lRISkgY6ciNgH7AO47oPrszvqQoKZ9CcNZiC48seSZs8vEHSs919ztwYZ6QwjIk6l+7PAz4GtwJnFoUa6P5s2PwVsaHv6+tRWO65p2LQaOjAkrZJ03eIysA14DjgI7Eqb7QIeT8sHgbvTtyW3Am+1DV26v8egfarw4HRNw6bRKEOSNcDP1TrlXgH8R0T8l6SjwGOS7gVeA76Qtj8E7ADmgXeBe/q9waSvwyj0fp3Dk/c8PLHmGjowIuJl4G+7tP8v8Nku7QHsGfb9cuWahk2TrK/0rAvXNGxaODBK4pqGTQMHRol8nYY1nQOjZB6eWJM5MMbAoWFNlX1gFD2gcjvwXNOwJso+MIp+DZnj15WuaVjTZB8YdefhiTWJA6MCDg1riuwDo641jE4ODWuC7AOjzjWMTg4Nq7vsA6NpHBpWZw6MCXBoWF01JjDqdmA5NKyOGhMYdahhdHJoWN00JjDqyqFhdeLAyIBDw+rCgZEJh4bVQfaBMU0Hh0PDcpd9YNSxmDkKh4blLPvAmEYODcuVAyNTDg3LUfaB0ZRfPhuGQ8Nyk31gNOmXz4bh0LCcZB8Y5tCwfDgwasKhYTnIPjCmuYbRyRML26RlHxjTXsPo5ImFbZKyDwy7mocnNikOjJpyaNgkZB8YrmH05pqGVS37wHANY3muaViVsg8M68/DE6tK38CQtF/SWUnPtbXdIOmwpBPpfnVql6SHJM1LOiZpc9tzdqXtT0jaNZ4fZ3o5NKwKRc4wfgRs72jbCxyJiI3AkfQY4HZgY7rtBh6GVsAADwC3AFuBBxZDph/XMIpzTcPGrW9gRMSTwLmO5p3AgbR8ALijrf2RaHkKuF7SWuA24HBEnIuIN4HDXB1C3d/fNYyBuKZh4zRsDWNNRJxOy28Aa9LyOuD1tu1OprZe7VeRtFvSnKS5CxfeGbJ7083DExuXkYueERFAaZ+6iNgXEVsiYsvKlavKetmp49CwcRg2MM6koQbp/mxqPwVsaNtufWrr1d6XaxjDc03DyjZsYBwEFr/p2AU83tZ+d/q25FbgrTR0eQLYJml1KnZuS219uYYxGtc0rEwr+m0g6SfAZ4AbJZ2k9W3Ht4HHJN0LvAZ8IW1+CNgBzAPvAvcARMQ5Sd8CjqbtvhkRnYVUG5OQYAa0EDADgS4PT2bPLxB0rHf4Wg99AyMi7uqx6rNdtg1gT4/X2Q/sH6h3VhqHhpUh+ys9XcMoj2saNqrsA8M1jHK5pmGjyD4wrHz+ytWG5cCYUg4NG4YDY4q5pmGDcmBMOdc0bBAODPPwxApzYBjg0LBiHBh2mUPD+nFg2BIODVuOA8Ou4tCwXhwY1pVDw7pxYFhPDg3r5MCwZTk0rJ0Dw/pyaNgiB4YV4tAwcGDYABwa5sCwgTg0ppsDwwbm0JheDgwbikNjOjkwbGgOjelTi8Do92Hzh3FyPAnPdKlFYPSb4NcTAE+WJ+GZHrUIDMufhyfTwYFhpXFoNF8tAsM1jPpwTaPZahEYrmHUi2sazVWLwLD68fCkmRwYNjYOjeapRWC4hlFfrmk0Sy0CwzWMenNNozlqERhWfx6eNIMDwyrj0Ki/voEhab+ks5Kea2v7hqRTkp5Jtx1t6+6XNC/pRUm3tbVvT23zkvYO0knXMJrDNY16K3KG8SNge5f270XEzel2CEDSJuBO4GPpOf8maVbSLPB94HZgE3BX2rYQ1zCaxTWN+uobGBHxJHCu4OvtBB6NiPci4hVgHtiabvMR8XJEnAceTdvalPLwpJ5GqWHcJ+lYGrKsTm3rgNfbtjmZ2nq1X0XSbklzkuYuXHhnhO5Z7hwa9TNsYDwMfBS4GTgNfKesDkXEvojYEhFbVq5cBbiG0WSuadTLimGeFBFnFpcl/QD4RXp4CtjQtun61MYy7f3fzzWMRgsJZkALATMQ6HJNY/b8AkHHev97T8xQZxiS1rY9/Dyw+A3KQeBOSddKugnYCPwWOApslHSTpGtoFUYPDt9taxoPT+qh7xmGpJ8AnwFulHQSeAD4jKSbgQBeBb4EEBHHJT0GPA9cBPZExKX0OvcBTwCzwP6IOF72D2P11vVM47zPNHLSNzAi4q4uzT9cZvsHgQe7tB8CDg3Uu0QRy344+q23+ugZGtfOMPueQ2PSanGlp2sY08XXaeSrFoFh08c1jTw5MCxbDo38NCIw/GFpLodGXhoRGK5hNJtDIx+NCAxrPodGHhwYVhsOjclzYFitODQmy4FhtePQmBwHhtWSQ2MyHBhWWw6N6jkwrNYcGtVyYFjtOTSq48CwRnBoVMOBYY3h0Bg/B4Y1ikNjvBwY1jgOjfFxYFgjOTTGw4FhjeXQKJ8DwxrNoVGu2gSG/zFtWP5jSeWpTWB4khwbhScWLkdtAsNsVB6ejM6BYVPFoTGa2gRGr388/6PaoFzTGF5tAqNXDcO1DRuGaxrDqU1gmJXNw5PBOTBsqjk0BlObwHANw8bFNY3iahMYrmHYOLmmUUxtAsNs3Dw86c+BYdbGobG8voEhaYOkX0t6XtJxSV9O7TdIOizpRLpfndol6SFJ85KOSdrc9lq70vYnJO0apKOuYVhVXNPorcgZxkXgqxGxCbgV2CNpE7AXOBIRG4Ej6THA7cDGdNsNPAytgAEeAG4BtgIPLIZMEa5hWJVc0+iub2BExOmI+H1afht4AVgH7AQOpM0OAHek5Z3AI9HyFHC9pLXAbcDhiDgXEW8Ch4HtZf4wZmXy8ORqA9UwJH0Y+DjwNLAmIk6nVW8Aa9LyOuD1tqedTG292s2y5dBYqnBgSPoA8FPgKxHx5/Z1ERFAKXtL0m5Jc5LmLlx450q7axg2Ia5pXFEoMCStpBUWP46In6XmM2moQbo/m9pPARvanr4+tfVqXyIi9kXElojYsnLlqivtrmHYBLmm0VLkWxIBPwReiIjvtq06CCx+07ELeLyt/e70bcmtwFtp6PIEsE3S6lTs3JbazGrBwxNYUWCbTwJfBJ6V9Exq+xrwbeAxSfcCrwFfSOsOATuAeeBd4B6AiDgn6VvA0bTdNyPiXBk/hFlVQoIZ0ELADAS6HBqz5xcIOtY37Ay4b2BExG+AXj/1Z7tsH8CeHq+1H9g/SAcXKaLrzu/VbjYuPUPj2hlm32t2aNTmSk/XMCwn01rTqE1gmOVmGmsaDgyzEUxbaNQ+MJrwj2D1Nk2hUfvAcA3DcjAtoVH7wDDLxTSEhgPDrERNDw0HhlnJmhwaDgyzMWhqaDgwzMakiaHhwDAbo6aFhgPDbMyaFBq1CYzcd6TZcpoSGrUJDF+gZXXXhNCoTWCYNUHdQ8OBYVaxOodGbQIjtx1nNoq6Tixcm8BwDcOapo6T8NQmMMyaqG7DEweG2YTVKTQcGGYZqEtNw4Fhlok61DQcGGYZyX144sAwy0zOoeHAMMtQrjUNB4ZZpnKsaTgwzDKW2/DEgWGWuZxCo1aB0bkzJn0Ri1lVcqlp1CowOn+fxL9fYtMkh5pGrQLDbNpNenjiwDCrmUmGRq0CwzUMs5ZJ1TT6BoakDZJ+Lel5ScclfTm1f0PSKUnPpNuOtufcL2le0ouSbmtr357a5iXtHbSzrmGYXTGJmsaKAttcBL4aEb+XdB3wO0mH07rvRcS/tm8saRNwJ/Ax4K+AX0n6m7T6+8DngJPAUUkHI+L5Mn4Qs2kUEsyAFgJmINDl4cns+QWCjvUj/ifbNzAi4jRwOi2/LekFYN0yT9kJPBoR7wGvSJoHtqZ18xHxMoCkR9O2DgyzEVQZGgPVMCR9GPg48HRquk/SMUn7Ja1ObeuA19uedjK19WrvfI/dkuYkzV248M7Sda5hmHVVVU2jcGBI+gDwU+ArEfFn4GHgo8DNtM5AvjN0L9pExL6I2BIRW1auXLV0nWsYZj1VUdMoUsNA0kpaYfHjiPgZQEScaVv/A+AX6eEpYEPb09enNpZpN7MSjHt4UuRbEgE/BF6IiO+2ta9t2+zzwHNp+SBwp6RrJd0EbAR+CxwFNkq6SdI1tAqjBwfqrZn1Nc7rNIqcYXwS+CLwrKRnUtvXgLsk3QwE8CrwJYCIOC7pMVrFzIvAnoi4BCDpPuAJYBbYHxHHB+msIpYkYudjM2spfKZxseTAiIjfAN2OykPLPOdB4MEu7YeWe15fATMLCyzMtk6MHBZmvRUJjVhZ8pAkJzEjQmLm0sKku2JWC/2GJ4OqVWCAQ8NsUGWGRu0CAxwaZoO6HBqLNYvF0HjfYBGgyPjiJ0lvAy9Ouh8F3Aj8adKdKMD9LFcd+lmkj38dER8q8mKFrsOYoBcjYsukO9GPpDn3szzuZ3nK7mMthyRmNhkODDMrLPfA2DfpDhTkfpbL/SxPqX3MuuhpZnnJ/QzDzDLiwDCzwrINjFHn/xxDf16V9Gyav3Qutd0g6bCkE+l+dWqXpIdS349J2jzGfu2XdFbSc21tA/dL0q60/QlJuyroY+VzwhboZ6/5a3Pbn5ObZzcisrvR+m3WPwIfAa4B/gBsmnCfXgVu7Gj7F2BvWt4L/HNa3gH8J60Lb28Fnh5jvz4NbAaeG7ZfwA3Ay+l+dVpePeY+fgP4py7bbkr/3tcCN6XPwWwVnwlgLbA5LV8HvJT6k9v+7NXPse/TXM8wtpLm/4yI88Di/J+52QkcSMsHgDva2h+JlqeA6zvmDylNRDwJnBuxX7cBhyPiXES8CRwGto+5j71cnhM2Il4BFueEHftnIiJOR8Tv0/LbwOL8tbntz1797KW0fZprYBSa/7NiAfxS0u8k7U5ta6I1STLAG8CatDzp/g/ar0n1t/Q5YcuipfPXZrs/VcE8u+1yDYwcfSoiNgO3A3skfbp9ZbTO/bL7jjrXfjGmOWHLoKvnr70sp/3ZpZ9j36e5BsZy84JOREScSvdngZ/TOp07szjUSPdn0+aT7v+g/aq8vxFxJiIuRcQC8AOu/CmKifZRXeavJcP92a2fVezTXAMjq/k/Ja1S6484IWkVsI3WHKYHgcUK+C7g8bR8ELg7VdFvBd5qO6WtwqD9egLYJml1Oo3dltrGRhnOCSt1n7+WzPZnr35Wsk/LrDKXeaNVgX6JVhX36xPuy0doVZD/ABxf7A/wF8AR4ATwK+CG1C5af+Xtj8CzwJYx9u0ntE4/L9Aag947TL+Af6BVDJsH7qmgj/+e+nAsfUjXtm3/9dTHF4Hbq/pMAJ+iNdw4BjyTbjsy3J+9+jn2fepLw82ssFyHJGaWIQeGmRXmwDCzwhwYZlaYA8PMCnNgmFlhDgwzK+z/AYDYFcJEIkQeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.imshow(F)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f397468e-a31c-42f1-a171-38b6e74a123a",
   "metadata": {},
   "source": [
    "# Gaussian quadrature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e0c919a8-e3ff-4a31-b62c-565b02c1f021",
   "metadata": {},
   "outputs": [],
   "source": [
    "gauss_weight = {1:np.array([2]),\n",
    "                2:np.array([1, 1]),\n",
    "                3:np.array([8/9, 5/9, 5/9]),\n",
    "                4:np.array([0.6521451548625462, 0.6521451548625462, 0.34785484513745385, 0.34785484513745385]),\n",
    "                5:np.array([0.5688888888888889, 0.47862867049936647, 0.47862867049936647, 0.23692688505618908, 0.23692688505618908])\n",
    "               }\n",
    "gauss_point = {1:np.array([0.0]),\n",
    "               2:np.array([0.5773502691896258, -0.5773502691896258]),\n",
    "               3:np.array([0.0, 0.7745966692414834, -0.7745966692414834]),\n",
    "               4:np.array([0.3399810435848563, -0.3399810435848563, 0.8611363115940526, -0.8611363115940526]),\n",
    "               5:np.array([0.0, 0.5384693101056831, -0.5384693101056831, 0.906179845938664, -0.906179845938664])\n",
    "              }\n",
    "\n",
    "def gauss_line_int(f, x1, y1, x2, y2, pt=2):\n",
    "    w = gauss_weight[pt]\n",
    "    x = gauss_point[pt]\n",
    "    \n",
    "    L = np.sqrt((x1-x2)**2+(y1-y2)**2)\n",
    "    ret = 0.5*L*np.sum(w*f(0.5*(x1-x2)*x+0.5*(x1+x2), 0.5*(y1-y2)*x+0.5*(y1+y2)))\n",
    "    return ret\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "220060c0-b92a-4c9f-bfe8-6ae8e9991f5a",
   "metadata": {},
   "source": [
    "# Boudary condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "4eca95e7-31db-4c88-ac97-2cf91c2c2ec3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.00017018929889045403\n",
      "0.0005106359745688018\n",
      "0.000851286911272871\n",
      "0.0011922783750701806\n",
      "0.0015337467655041684\n",
      "0.001875828675531734\n",
      "0.0022186609414643547\n",
      "0.0025623807010893285\n",
      "0.0029071254483526315\n",
      "0.0032530330817471718\n",
      "0.0036002419756281778\n",
      "0.0039488910155665745\n",
      "0.004299119668485705\n",
      "0.0046510680204725425\n",
      "0.005004876861422321\n",
      "0.005360687709588138\n",
      "0.0057186429353101335\n",
      "0.006078885648835355\n",
      "0.006441560047168222\n",
      "0.006806811118592225\n",
      "0.00717478499442384\n",
      "0.007545628930895713\n",
      "0.00791949120386374\n",
      "0.00829652134666607\n",
      "0.008676870237632282\n",
      "0.009060689904884294\n",
      "0.00944813405974647\n",
      "0.009839357604872895\n",
      "0.010234516965008355\n",
      "0.010633770272469626\n",
      "0.011037277300349074\n",
      "0.011445199314399577\n",
      "0.011857699524095584\n",
      "0.012274943113326045\n",
      "0.012697096771614718\n",
      "0.013124329388561392\n",
      "0.013556812030882824\n",
      "0.01399471743311695\n",
      "0.014438221018964803\n",
      "0.014887500103169648\n",
      "0.015342734270031117\n",
      "0.01580410591451538\n",
      "0.01627179926502505\n",
      "0.016746001570636077\n",
      "0.017226902501254242\n",
      "0.017714694485426234\n",
      "0.018209572536610173\n",
      "0.018711734717347023\n",
      "0.01922138175927826\n",
      "0.01973871752926601\n",
      "0.01987025253150939\n",
      "0.019605476118493104\n",
      "0.019332857843407558\n",
      "0.019052506726703597\n",
      "0.018764534846811325\n",
      "0.018469057470518777\n",
      "0.018166192588951576\n",
      "0.01785606158587657\n",
      "0.017538788359261072\n",
      "0.01721449981733763\n",
      "0.01688332571719309\n",
      "0.016545398518871372\n",
      "0.016200853340883236\n",
      "0.015849828062089462\n",
      "0.015492463156321159\n",
      "0.015128901288309689\n",
      "0.014759288201310292\n",
      "0.014383771587278788\n",
      "0.014002501598068033\n",
      "0.013615630801575162\n",
      "0.0132233138879297\n",
      "0.01282570791818453\n",
      "0.012422971846871063\n",
      "0.012015266743887294\n",
      "0.011602755667215062\n",
      "0.011185603587395131\n",
      "0.010763977446678511\n",
      "0.010338045865478119\n",
      "0.009907979231348132\n",
      "0.009473949533862571\n",
      "0.009036130397640797\n",
      "0.008594696894258536\n",
      "0.00814982566711406\n",
      "0.007701694579221939\n",
      "0.0072504829475318385\n",
      "0.006796371175597961\n",
      "0.006339540978579191\n",
      "0.005880175030376383\n",
      "0.005418457145734075\n",
      "0.00495457188627182\n",
      "0.004488704885000259\n",
      "0.004021042494136375\n",
      "0.003551771682359216\n",
      "0.0030810802400987676\n",
      "0.0026091564139573502\n",
      "0.002136188948409761\n",
      "0.0016623670292453356\n",
      "0.0011878802064373163\n",
      "0.000712918233588636\n",
      "0.000237671098376211\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "# def dir1(x, y):\n",
    "#     return x*(1-x)\n",
    "\n",
    "# def dir2(x, y):\n",
    "#     return -y*(1-y)\n",
    "\n",
    "# def dir3(x, y):\n",
    "#     return x*(1-x)\n",
    "\n",
    "# def dir4(x, y):\n",
    "#     return -y*(1-y)\n",
    "\n",
    "def dir1(x, y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def dir2(x, y):\n",
    "    return np.sinh(y)/np.sinh(1)\n",
    "\n",
    "def dir3(x, y):\n",
    "    return np.sin(x)/np.sin(1)\n",
    "\n",
    "def dir4(x, y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "\n",
    "du = np.zeros(N)\n",
    "\n",
    "for edge in dir_edge_1:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir1, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L\n",
    "    print(up)\n",
    "\n",
    "for edge in dir_edge_2:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir2, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L\n",
    "    print(up)\n",
    "\n",
    "for edge in dir_edge_3:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir3, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L\n",
    "    print(up)\n",
    "\n",
    "for edge in dir_edge_4:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir4, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L\n",
    "    print(up)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "de68892b-9135-471f-8dc5-2cdef4f22f57",
   "metadata": {},
   "outputs": [],
   "source": [
    "# with np.printoptions(precision=2, suppress=True):\n",
    "#     print(F)\n",
    "#     print(F[0:-1,0:-1])\n",
    "#     print(du[-1])\n",
    "#     print(du)\n",
    "#     print(s)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "91bcc7e5-0b43-4991-a0ab-86921cfbc58b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# s = np.zeros(N)\n",
    "\n",
    "# s[0:-1] = np.linalg.inv(F[0:-1,0:-1])@du[0:-1]\n",
    "\n",
    "\n",
    "FF = F[:]\n",
    "ddu = du[:]\n",
    "FF[-1,:] = 0\n",
    "FF[:,-1] = 0\n",
    "FF[-1,-1] = 1\n",
    "ddu[-1] = 0\n",
    "\n",
    "s = np.linalg.inv(FF) @ ddu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "39c9e64e-3169-42cd-be4c-afc135fa7f77",
   "metadata": {},
   "outputs": [],
   "source": [
    "E = len(element)\n",
    "q_h = np.zeros((E,2)) \n",
    "\n",
    "for i, e in enumerate(element):\n",
    "    x1,y1 = node[e][0]\n",
    "    x2,y2 = node[e][1]\n",
    "    x3,y3 = node[e][2]\n",
    "    \n",
    "    A_e = 1/2*np.linalg.det(np.c_[node[e], np.ones((3,1))])\n",
    "    \n",
    "    Phi_e = 1/(2*A_e)*np.array([[x2-x3,x3-x1,x1-x2],\n",
    "                                [y2-y3,y3-y1,y1-y2]])\n",
    "    \n",
    "    q_e = Phi_e @ s[e]\n",
    "    q_h[i] = q_e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "61b4bd15-257e-4be2-ae82-3c06837546a2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<visvis.wibjects.colorWibjects.Colorbar object at 0x00000205439122B0>\n",
      "ERROR calling '_OnAxesPositionChange':\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\wibjects\\colorWibjects.py\", line 622, in _OnAxesPositionChange\n",
      "    self.position.x = axes.position.width+5\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\misc.py\", line 217, in fsetWithDraw\n",
      "    fset(self, *args)\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 1126, in fset\n",
      "    self._Update()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 928, in _Update\n",
      "    self._CalculateInPixels()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 997, in _CalculateInPixels\n",
      "    raise Exception(\"Can only calculate the position in pixels\"+\n",
      "Exception: Can only calculate the position in pixels if the owner has a parent!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "cfv.figure()\n",
    "\n",
    "# Draw the mesh.\n",
    "\n",
    "cfv.drawElementValues(\n",
    "    ev=q_h[:,0],\n",
    "    coords=coords,\n",
    "    edof=edof,\n",
    "    \n",
    "    dofs_per_node=mesh.dofsPerNode,\n",
    "    el_type=mesh.elType,\n",
    "    title=\"Example 01\"\n",
    "        )\n",
    "cfv.showAndWait()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "7aebb1dc-ad69-4ce2-84fc-81008e951af4",
   "metadata": {},
   "outputs": [],
   "source": [
    "f = lambda x:x**2/2-x**3/3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "8e1a62a6-d94b-4cd3-8a39-26240c39f14a",
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = np.meshgrid(np.arange(0,1.01,0.01), np.arange(0,1.01,0.01))\n",
    "z = np.cos(x)*np.sinh(y)/(np.sin(1)*np.sinh(1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "21addc22-ed08-431e-a0c2-a4528c1c4e5a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAFpCAYAAABu98hvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAVe0lEQVR4nO3ca4yldX3A8e+PpUCwKzSuRsOulbarlWATLUGNSbXxkpUX7AtbC4S0tNSNtpimtk1obKjBN9pGG02ourXUS6J4aWImcZXeJCRG6G6iVdlUs6DAoBYviEktrKu/vjjnwNlhZs4z5zzX//P9JJOcy39mfvvfme8++5xLZCaSpGE7o+sBJEmrM+aSVABjLkkFMOaSVABjLkkFMOaSVICFMY+ImyPiwYj46hb3R0S8OyJORMSXI+IF9Y8pSdpOlSPzDwAHtrn/1cD+6cch4D2rjyVJ2omFMc/M24EfbLPkIPChnLgDOD8inlHXgJKkxeo4Z34BcP/c9fXpbZKklpzZ5jeLiENMTsVw7rnx6xf+cqvfXgLgoZ+e2/UIvfejU+d0PUIrTp7sV4NOfvOB72XmU5f53Dr+JA8A++au753e9gSZeRg4DHDxr52V//zpPTV8e2mxT/7Ix+Wr+Nfv/GrXI7Tmm+tLNbNR915z/b3Lfm4dMV8DrouIW4AXAg9n5rdr+LrSSgz4zhjyYVsY84j4KPAyYE9ErAN/DfwcQGa+FzgCXAacAH4M/H5Tw0pVGPGdMeJlWBjzzLxywf0J/HFtE0lLMODLMeTl6NfZf2kHDPjyxhRxKD/kYMw1QEZ8NYa8TMZcg2DAVze2iMN4Qg7GXD1mwOsxxojDuEIOxlw9Y8DrM9aIw/hCDsZcHTPezRhryMcY8RljrtYZ8OaMNeIw7pCDMVdLDHizxhxxMORgzNUgA94OQ27IwZirRsa7XWOPOBjyecZcSzPe3TDiRnwzxlyVGe9uGfEJQ745Y64tGe9+MOKPM+RbM+Z6jPHuFyN+utJDfvZ9Z630+cZ8xIx3Pxnx05UecVg95GDMR8Nw958RfyJDXp0xL5TxHg4jvjlDvjPGfOCM9nAZ8a2NIeR1M+YDYrjLYMS3NqaI13lUDsa8l4x2mYz49gz5aox5Rwz2eBjxxQz56ox5gwz2uBnxxcYUcWgu5GDMl2aotRUjXo0hr5cx38BIa1lGvDpDXr8iY26Q1RYDvjNjizi0E3LoMOYP/fRco6vBMuI7Z8ibVeSRudQUI75zY4w4tBtyMObSQgZ8eYa8PcZc2oIRX95YIw7dhByMuXQaA746Q94NYy5hxOsw5ohDtyEHY64RM+D1MeTdhhyMuUbGgNdr7BGHfoQcjLlGwojXz5D3J+RgzFUwA94MIz7Rp5CDMVdhDHizDPlE30IOxlwFMODNM+KP62PIwZhroAx4O4z46foacjDmGhAD3i5Dfro+hxyMuXrOgLfPiD9R30MOxlw9ZMC7YcSfaAgRnzHm6pzx7p4hf6IhhRyMuTpiwPvBiG9uaCEHY66WGO9+MeJbG2LIwZirQQa8f4z49oYacjDmqpHx7jdDvr0hhxyMuVZgvIfBiC829JCDMdcOGO9hMeKLlRDxGWOuLRnvYTLi1ZQUcjDmmjLcw2fEqyst5GDMR8lwl8WI70yJIQdjXjzDXS4jvnOlhhyMeVEM9zgY8Z0rOeIzxnyAjPY4GfHljCHkYMx7zWgLjPgqxhJyMOadM9jaihFfzZhCDhVjHhEHgHcBu4D3Z+bbNtz/TOCDwPnTNddn5pF6Rx0mY62dMuKrGVvEZxbGPCJ2ATcBrwTWgaMRsZaZx+eW/RXw8cx8T0RcBBwBntXAvL1iqFUnI766sYYcqh2ZXwqcyMx7ACLiFuAgMB/zBJ48vXwe8K06h2yLcVYXjHg9xhxyqBbzC4D7566vAy/csOYtwL9ExBuBJwGvqGW6HTDEGhIDXp+xR3ymrgdArwQ+kJnviIgXAx+OiIsz82fziyLiEHAI4Oyn7TbAGh0jXi9D/rgqMX8A2Dd3fe/0tnnXAgcAMvMLEXEOsAd4cH5RZh4GDgPsfs7Tc8mZpcEx4vUz5KerEvOjwP6IuJBJxK8Artqw5j7g5cAHIuK5wDnAd+scVBoaA94MI765hTHPzFMRcR1wK5OnHd6cmXdFxI3AscxcA/4M+IeI+FMmD4Zek5keeWuUjHhzDPnWKp0znz5n/MiG226Yu3wceEm9o0nDYsSbY8QX8xWg0goMePMMeTXGXFqCEW+HIa/OmEsVGfD2GPGdM+bSNgx4+wz5coy5tAkj3j4jvhpjLk0Z8O4Y8tUZc42aAe+WEa+PMdfoGPB+MOT1MuYaBQPeH0a8GcZcxTLg/WLEm2XMVRQD3k+GvHnGXINnwPvLiLfHmGtwjPcwGPKd2X3vam80a8w1CAZ8OIz4zq0acjDm6injPTxGfDl1hByMuXrEgA+TEV9eXSEHY64OGe/hM+TLqzPkYMzVIuNdDiO+mrpDDsZcDTLe5THiq2si5GDMVSPjXS4jXo+mQg7GXEsy3ONhyOvRZMjBmKsi4z0+Rrw+TYccjLk2YbjHzYjXq42QgzEfPcOtGSNev7ZCDsZ8NIy2tmLE69dmxGeMeWGMtqoy4s3oIuRgzAfLaGtZRrw5XYUcjHmvGWzVyYg3q8uQgzHvlLFWG4x487oOORjzxhhqdc2It6MPIQdjvmNGWn1mwNvTl4jPjDLmBlmlMeLt6lvIocOYnzx5plGVVmTE29fHkMNIj8yloTPi3ehryMGYS4NixLvR54jPGHOp5wx4t4YQcjDmUm8Z8e4NJeRgzKVeMeD9MKSIzxhzqQeMeH8MMeRgzKXOGPD+GWrIwZhLrTPi/TPkiM8Yc6kFBry/Sgg5GHOpMQa830qJ+Iwxl2pmxPuvtJCDMZdqYcCHo8SQgzGXlmbAh6XUiM8Yc2kHDPgwlR5yMObSQgZ8uMYQ8RljLm3CgA/fmEIOxlx6jAEvw9giPmPMNVrGuyxjjfiMMdeoGPAyjT3kYMxVOONdNiP+OGOu4hjwcTDkpzPmGjzjPS5GfHPGXINjvMfJiG/PmKv3jLcM+WJnVFkUEQci4msRcSIirt9izWsj4nhE3BURH6l3TI3J2fedddqHxmv3vWnIK1p4ZB4Ru4CbgFcC68DRiFjLzONza/YDfwm8JDMfioinNTWwymOwtZEB37kqp1kuBU5k5j0AEXELcBA4PrfmdcBNmfkQQGY+WPegKofx1laM+PKqxPwC4P656+vACzeseTZARHwe2AW8JTM/u/ELRcQh4BDArqecv8S4GhrDraoM+WrqegD0TGA/8DJgL3B7RDwvM384vygzDwOHAc6+cK9/c4Ux3FqGEa9HlZg/AOybu753etu8deDOzPwJ8I2I+DqTuB+tZUr1juHWqox4varE/CiwPyIuZBLxK4CrNqz5FHAl8E8RsYfJaZd7apxTHTHaqpsRb8bCmGfmqYi4DriVyfnwmzPzroi4ETiWmWvT+14VEceBnwJ/kZnfb3Jw1c9wq0lGvFmVzpln5hHgyIbbbpi7nMCbph/qMYOtthnxdvgK0EIZbXXNiLfLmA+UsVZfGfFuGPOeMtYaGiPeLWPeAUOtkhjxfjDmNTHQGhsj3i/GfAOjLG3PiPfToGJuaKVuGPD+6yzmcTKMs9RzRnw4BnVkLqkdRnx4jLmkxxjx4TLm0sgZ8DIYc2mkjHhZjLk0Iga8XMZcGgEjXj5jLhXKgI+LMZcKY8THyZhLBTDgMubSQBlwzTPm0oAYcG3FmEsDYMS1iDGXesqAj8t5dz+60ucbc6lHDPj4rBrxGWMudcyAj1ddIQdjLrXOeAvqDTkYc6kVBlwzdUd8xphLDTHg2qipkIMxl2pjvLWVJiM+Y8ylFRhwLdJGyMGYSztivFVVWxGfMebSNoy3ltF2yMGYS09gwLWsLiI+Y8w1esZbdegy5GDMNTKGW3XrOuIzxlxFM95qSl8iPmPMVRTjrTb0LeRgzDVghltt62PEZ4y5BsFwq0t9jviMMVfvGG71yRBCDsZcHTPc6quhRHzGmKsVRltDMbSIzxhz1cpoa6iGGvEZY66lGG2VYugRnzHm2pLBVslKifiMMR8xY60xKi3iM8a8YMZaelypEZ8x5gNlqKVqSo/4jDHvGSMt1WMsEZ8x5g0yzFL7xhbxGWO+DWMsDcdYIz7T+5gbVEnbGXvEZzqL+a6ThlrS8oz46Xp/ZC5JMwZ8a8ZcUu8Z8cWMuaTeMuLVGXNJvWLAl2PMJfWCEV/NGVUWRcSBiPhaRJyIiOu3WfeaiMiIuKS+ESWV6ry7H33sQ6tZeGQeEbuAm4BXAuvA0YhYy8zjG9btBv4EuLOJQSWVw3jXr8qR+aXAicy8JzNPArcABzdZ91bg7cAjNc4nqRAehTerSswvAO6fu74+ve0xEfECYF9mfnq7LxQRhyLiWEQcO/V//7vjYSUNiwFvz8oPgEbEGcA7gWsWrc3Mw8BhgCc9dZ8v/5QKZbzbVyXmDwD75q7vnd42sxu4GLgtIgCeDqxFxOWZeayuQSX1mwHvVpWYHwX2R8SFTCJ+BXDV7M7MfBjYM7seEbcBf27IpfIZ8P5YGPPMPBUR1wG3AruAmzPzroi4ETiWmWtNDympPwx4P1U6Z56ZR4AjG267YYu1L1t9LEl9YsD7z1eAStqUAR8WYy7pMQZ8uIy5NGLGuxzGXBoZA14mYy4VzniPgzGXCmTAx8eYSwUw3jLm0gAZb21kzKUBMN5axJhLPWS8tVPGXOoB461VGXOpZYZbTTDmUsOMt9pgzKUaGW51xZhLSzLc6hNjLlVguNV3xlyaY7Q1VMZco2W4VRJjruIZbY2BMVcRDLbGzphrUIy2tDljrl4x1tJyjLlaZ7Cl+hlz1cpQS90w5qrMUEv9ZcxlpKUCGPMCGWdpfIx5zxlmSVUY8wYZYkltGV3MDaykEnUW812PpGGVpJqc0fUAkqTVGXNJKoAxl6QCGHNJKoAxl6QCGHNJKoAxl6QeOOu/11f6/NG9aEiS+mLVgM8z5pLUsjojPmPMJakFTQR8njGXpAY1HfEZYy5JNWsr4POMuSTVoIuAzzPmkrSkrgM+z5hL0g70KeDzjLkkLdDXgM8z5pK0wRDivZExlySGGfB5xlzSKA093hsZc0mjUVrA5xlzScUqOd4bGXNJxRhTvDcy5pIGa8zx3siYSxoM4701Yy6plwz3zlSKeUQcAN4F7ALen5lv23D/m4A/BE4B3wX+IDPvrXlWSQUz3qtZGPOI2AXcBLwSWAeORsRaZh6fW/ZF4JLM/HFEvAH4G+B3mhhY0vAZ7vpVOTK/FDiRmfcARMQtwEHgsZhn5ufm1t8BXF3nkJKGy3C3o0rMLwDun7u+Drxwm/XXAp9ZZShJw2S4u1PrA6ARcTVwCfDSLe4/BBwCOPvs8+v81pJaZLT7p0rMHwD2zV3fO73tNBHxCuDNwEsz89HNvlBmHgYOAzx5997c8bSSWmW0h6NKzI8C+yPiQiYRvwK4an5BRDwfeB9wIDMfrH1KSY0y2sO3MOaZeSoirgNuZfLUxJsz866IuBE4lplrwN8CPw98IiIA7svMyxucW9ISjHa5Kp0zz8wjwJENt90wd/kVNc8laUkGe5x8Bag0QAZbGxlzqYeMtXbKmEsdMNaqmzGXamao1QVjLlVkpNVnxlyjZ6RVAmOu4hhnjZExV68ZZqkaY67GGWSpecZcWzLC0nAY8wIYXUnGvAbGVFLXOot5PHLSCEpSTc7oegBJ0uqMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgEqxTwiDkTE1yLiRERcv8n9Z0fEx6b33xkRz6p9UknSlhbGPCJ2ATcBrwYuAq6MiIs2LLsWeCgzfwX4O+DtdQ8qSdpalSPzS4ETmXlPZp4EbgEOblhzEPjg9PIngZdHRNQ3piRpO1VifgFw/9z19eltm67JzFPAw8BT6hhQkrTYmW1+s4g4BByaXn30s9/5+6+2+f17aA/wva6H6AH3wT0A9wDgOct+YpWYPwDsm7u+d3rbZmvWI+JM4Dzg+xu/UGYeBg4DRMSxzLxkmaFL4R5MuA/uAbgHMNmDZT+3ymmWo8D+iLgwIs4CrgDWNqxZA35vevm3gP/IzFx2KEnSziw8Ms/MUxFxHXArsAu4OTPviogbgWOZuQb8I/DhiDgB/IBJ8CVJLal0zjwzjwBHNtx2w9zlR4Df3uH3PrzD9SVyDybcB/cA3ANYYQ/CsyGSNHy+nF+SCtB4zH0rgEp78KaIOB4RX46If4+IX+xiziYt2oO5da+JiIyI4p7VUGUPIuK105+FuyLiI23P2IYKvw/PjIjPRcQXp78Tl3UxZ1Mi4uaIeDAiNn1qdky8e7o/X46IF1T6wpnZ2AeTB0zvBn4JOAv4L+CiDWv+CHjv9PIVwMeanKntj4p78JvAudPLbxjjHkzX7QZuB+4ALul67g5+DvYDXwR+YXr9aV3P3dE+HAbeML18EfDNrueueQ9+A3gB8NUt7r8M+AwQwIuAO6t83aaPzH0rgAp7kJmfy8wfT6/eweS5/CWp8nMA8FYm7+vzSJvDtaTKHrwOuCkzHwLIzAdbnrENVfYhgSdPL58HfKvF+RqXmbczedbfVg4CH8qJO4DzI+IZi75u0zH3rQCq7cG8a5n8q1yShXsw/a/kvsz8dJuDtajKz8GzgWdHxOcj4o6IONDadO2psg9vAa6OiHUmz6J7Yzuj9cZOmwG0/HJ+bS8irgYuAV7a9SxtiogzgHcC13Q8StfOZHKq5WVM/nd2e0Q8LzN/2OVQHbgS+EBmviMiXszkNSwXZ+bPuh6sz5o+Mt/JWwGw3VsBDFiVPSAiXgG8Gbg8Mx9taba2LNqD3cDFwG0R8U0m5wnXCnsQtMrPwTqwlpk/ycxvAF9nEveSVNmHa4GPA2TmF4BzmLxvy1hUasZGTcfctwKosAcR8XzgfUxCXuJ50m33IDMfzsw9mfmszHwWk8cNLs/Mpd+nooeq/C58islRORGxh8lpl3tanLENVfbhPuDlABHxXCYx/26rU3ZrDfjd6bNaXgQ8nJnfXvhZLTxyexmTI4y7gTdPb7uRyS8rTP6iPgGcAP4T+KWuH23uYA/+Dfgf4EvTj7WuZ257DzasvY3Cns1S8ecgmJxuOg58Bbii65k72oeLgM8zeabLl4BXdT1zzX/+jwLfBn7C5H9j1wKvB14/93Nw03R/vlL1d8FXgEpSAXwFqCQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgGMuSQVwJhLUgH+H1bfS8wmT6VSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(6,6))\n",
    "plt.contourf(x, y, z)\n",
    "plt.xlim([0,1])\n",
    "plt.ylim([0,1])\n",
    "plt.gca().set_aspect('equal', adjustable='box')\n",
    "# plt.axis('equal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "dd938826-63d8-4824-8c0f-fe4c315ace2d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[   1,    5,  200],\n",
       "       [ 200,    5,  201],\n",
       "       [ 200,  201,  199],\n",
       "       ...,\n",
       "       [2601,  101,  102],\n",
       "       [2601,  102,  103],\n",
       "       [ 103,  102,    3]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "edof"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e897cb5-25ff-4485-9a00-07431934432a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "EFEM",
   "language": "python",
   "name": "efem"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
